Violations of physical and psychological expectations in the human adult brain - Experiment 2 Univariate analysis
Shari Liu
2023-08-02
1 Overview
This script visualizes and analyzes the mean amplitude of response
(meanbeta) to expected and unexpected events from three
“domains” - psychology-action, psychology-environment, and physics. Our
pre-registered confirmatory analyses center around the psychology-action
and physics videos.
- Part 1a: Domain and event effects 8 focal regions. This is the univariate confirmatory analysis for the project.
- Part 1b: Checking robustness of results depending on ROI size (pre-registered as 100 voxels)
- Part 2: Repeating the univariate analysis in a larger set of domain-general and domain-specific regions, with more conservative significance thresholds given the # of regions we’re exploring per category.
- Part 3: Exploratory univariate analysis on VOE effect in
action-environment scenarios.
- Part 4: Similar to confirmatory analysis, except that we use the VOE events from runs 3-4 to select ROIs. Function of this analysis is to check for whether the confusing STS results we obtained is due to conceptual mismatch between the DOTS localizer and the psych-action events.
- Part 5: We model betas for each scenario and event, rather than folding across scenarios / tasks. This allows us to test for effects of visual statistics (Part 5a), which were computed over each video individually. Also allows us to do per-task analyses (Part 5b).
- Part 6: Compute and plot Dice coefficient between two focal regions (bilateral APC and L/RSMG).
- Part 7: Region by region univariate effect size (analogue to confirmatory MVPA analysis that is rendered uninterpretable, given the lack of reliable VOE spatial patterns across runs)
knitr::opts_chunk$set(echo = TRUE)
options(scipen = 1, digits = 3)
library(pacman)
pacman::p_load(tidyverse,
patchwork,
papaja,
ggforce,
conflicted,
here,
cowplot,
lme4,
lmerTest,
effects,
effectsize,
EMAtools,
performance,
parameters,
BayesFactor,
gt,
insight,
sjPlot,
lsmeans,
ggrepel,
GGally,
nptest,
boot,
boot.pval,
car)
conflict_prefer("select", "dplyr")## [conflicted] Will prefer dplyr::select over any other packageconflict_prefer("filter", "dplyr")## [conflicted] Will prefer dplyr::filter over any other packageconflict_prefer("rename", "dplyr")## [conflicted] Will prefer dplyr::rename over any other packageconflict_prefer("mutate", "dplyr")## [conflicted] Will prefer dplyr::mutate over any other packageconflict_prefer("summarise", "dplyr")## [conflicted] Will prefer dplyr::summarise over any other packageconflict_prefer("lmer", "lmerTest")## [conflicted] Will prefer lmerTest::lmer over any other packageconflict_prefer("here", "here")## [conflicted] Will prefer here::here over any other packagesource(here("helper_funs.R"))
options(contrasts = c("contr.sum", "contr.poly"))
theme_set(theme_cowplot(10))
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] car_3.0-8 boot.pval_0.4 boot_1.3-25
## [4] nptest_1.0-3 GGally_2.1.2 ggrepel_0.8.2
## [7] lsmeans_2.30-0 emmeans_1.4.8 sjPlot_2.8.10
## [10] insight_0.15.0 gt_0.8.0 BayesFactor_0.9.12-4.4
## [13] coda_0.19-3 parameters_0.16.0 performance_0.8.0
## [16] EMAtools_0.1.3 effectsize_0.6.0.1 effects_4.1-4
## [19] carData_3.0-4 lmerTest_3.1-2 lme4_1.1-23
## [22] Matrix_1.2-18 cowplot_1.0.0 here_1.0.1
## [25] conflicted_1.1.0 ggforce_0.3.4 papaja_0.1.0.9997
## [28] tinylabels_0.2.3 patchwork_1.1.2 forcats_0.5.0
## [31] stringr_1.5.0 dplyr_1.0.9 purrr_0.3.4
## [34] readr_1.3.1 tidyr_1.2.0 tibble_3.1.7
## [37] ggplot2_3.4.1 tidyverse_1.3.0 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.4.1 plyr_1.8.6
## [4] splines_4.0.2 TH.data_1.0-10 digest_0.6.31
## [7] htmltools_0.5.4 fansi_1.0.3 magrittr_2.0.3
## [10] memoise_2.0.1 openxlsx_4.1.5 modelr_0.1.8
## [13] sandwich_2.5-1 colorspace_2.0-3 blob_1.2.1
## [16] rvest_0.3.6 mitools_2.4 rbibutils_2.2.8
## [19] haven_2.4.3 xfun_0.36 crayon_1.5.1
## [22] jsonlite_1.8.4 survival_3.1-12 zoo_1.8-8
## [25] glue_1.6.2 polyclip_1.10-0 gtable_0.3.0
## [28] MatrixModels_0.4-1 sjstats_0.18.0 sjmisc_2.8.5
## [31] abind_1.4-5 scales_1.2.0 mvtnorm_1.1-1
## [34] DBI_1.1.2 ggeffects_0.16.0 Rcpp_1.0.8
## [37] xtable_1.8-4 foreign_0.8-80 survey_4.0
## [40] datawizard_0.2.3 httr_1.4.2 RColorBrewer_1.1-3
## [43] ellipsis_0.3.2 pkgconfig_2.0.3 reshape_0.8.8
## [46] farver_2.1.0 nnet_7.3-14 sass_0.4.4
## [49] dbplyr_1.4.4 utf8_1.2.2 tidyselect_1.1.2
## [52] rlang_1.0.6 munsell_0.5.0 cellranger_1.1.0
## [55] tools_4.0.2 cachem_1.0.6 cli_3.5.0
## [58] generics_0.1.2 sjlabelled_1.1.7 broom_0.8.0
## [61] evaluate_0.19 fastmap_1.1.0 yaml_2.3.6
## [64] knitr_1.41 fs_1.5.2 zip_2.2.0
## [67] pbapply_1.6-0 nlme_3.1-148 xml2_1.3.2
## [70] compiler_4.0.2 rstudioapi_0.13 curl_4.3
## [73] reprex_0.3.0 statmod_1.4.34 tweenr_1.0.2
## [76] bslib_0.4.2 stringi_1.7.8 lattice_0.20-41
## [79] nloptr_1.2.2.2 vctrs_0.5.1 pillar_1.7.0
## [82] lifecycle_1.0.3 Rdpack_2.3.1 jquerylib_0.1.4
## [85] estimability_1.3 data.table_1.13.0 R6_2.5.1
## [88] rio_0.5.16 codetools_0.2-16 MASS_7.3-51.6
## [91] assertthat_0.2.1 rprojroot_2.0.3 withr_2.5.0
## [94] multcomp_1.4-13 bayestestR_0.11.5 hms_0.5.3
## [97] DataCombine_0.2.21 grid_4.0.2 minqa_1.2.4
## [100] rmarkdown_2.19.2 numDeriv_2016.8-1.1 lubridate_1.7.9regioninfo <- read.csv(here("input_data/manyregions_info.csv"))
DT::datatable(regioninfo %>% dplyr::arrange(ROI_category, focal_region), options = list(scrollX = TRUE))domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")univariate_data <-
read_rds(here("outputs/univariate_data/study2_univariate_data_allregions_alltopN.Rds"))
univariate_data$event <- relevel(univariate_data$event, ref = "unexp") # for some reason this is required for the model coefficients to reflect the intuitive def that unexpected > expected is a positive coefficient
univariate_summary_domain <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain.Rds"))
univariate_summary_domain_runs12 <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_runs12.Rds"))
univariate_summary_domain_allruns <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_allruns.Rds"))2 PART 1a: Focal regions - confirmatory analysis
focal_data <- univariate_data %>%
filter(focal_region == 1)
focal_summary_domain <- univariate_summary_domain %>%
filter(focal_region == 1)
focal_summary_domain_runs12 <- univariate_summary_domain_runs12 %>%
filter(focal_region == 1)
# focal_summary_task_runs12 <- univariate_summary_task_runs12 %>%
# filter(focal_region == 1)
regions <- levels(as.factor(focal_data$ROI_name_final))2.1 Plots
First, here are plots of betas, run by run, to physics and psychology events (familiarization, expected, unexpected).
plot_univar_runbyrun("superiortemporal_L", -1, 4) | plot_univar_runbyrun("superiortemporal_R", -1, 4)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("supramarginal_L", -1, 4) | plot_univar_runbyrun("supramarginal_R", -1, 4)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("antParietal_bilateral", -1, 6) | plot_univar_runbyrun("precentral_inferiorfrontal_R", -1, 6)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("V1_bilateral", -1, 11) | plot_univar_runbyrun("MT_bilateral", -1, 11)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
2.1.1 By domain (input from domain / event-level betas)
These plots go into Figure 3 of the main text.
# DS_runs12_plots <- plot_univar_runs12("superiortemporal", 0, 4) |
# plot_univar_runs12("supramarginal", 0, 4)
#
# DG_runs12_plots <- plot_univar_runs12("MD", 0, 6) | plot_univar_runs12("early_visual", 0, 11)
#
# DS_runs12_plots / DG_runs12_plots
(plot_univar_runs12("superiortemporal_L", 0, 4) | plot_univar_runs12("superiortemporal_R", 0, 4) )## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
(plot_univar_runs12("supramarginal_L", 0, 4) | plot_univar_runs12("supramarginal_R", 0, 4))## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
(plot_univar_runs12("antParietal_bilateral", 0, 6) | plot_univar_runs12("precentral_inferiorfrontal_R", 0, 6))## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_runs12("V1_bilateral", 0, 11) | plot_univar_runs12("MT_bilateral", 0, 11)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
2.1.2 Action-env (input from domain / event-level betas)
summary_domain_allruns <- read_rds(here("outputs/descriptive_summaries/study2_univariate_summary_domain_allruns.Rds"))Univariate response to psych_env events, over all runs, for all focal regions.
STS <- plot_both_allruns("superiortemporal_", -1, 10)
SMG <- plot_both_allruns("supramarginal_", -1, 10)
V1 <- plot_both_allruns("V1_bilateral", -1, 10)
MT <- plot_both_allruns("MT_bilateral", -1, 10)
APC <- plot_both_allruns("antParietal_bilateral", -1, 10)
RFC <- plot_both_allruns("precentral_inferiorfrontal_R", -1, 10)
STS | SMG | V1 | MT | APC | RFC + plot_layout(widths = c(2, 2, 1, 1, 1, 1))## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.2.1.3 By task (input from single betas)
study2_univariate_summary_task_singlebetas <- readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_task_singlebetas.Rds")) %>%
filter(focal_region == "1")regions <- levels(as.factor(study2_univariate_summary_task_singlebetas$ROI_name_final))
plot_univar <- function(region, ymin, ymax) {
plotobject <- ggplot(data = study2_univariate_summary_task_singlebetas %>%
filter(str_detect(ROI_name_final, region)),
aes(x = event, y = meanbeta, fill = task)) +
geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
geom_errorbar(
aes(ymin = meanbeta - se, ymax = meanbeta + se),
position = position_dodge(width = .9),
width = .2,
colour = "black"
) +
theme_cowplot(10) +
facet_wrap( ~ ROI_name_final + factor(task, c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint")) + extracted_run_number, nrow = 2) +
# scale_fill_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
ylab("Average beta (amplitude)") +
xlab("Event") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
),
legend.position = "none") +
coord_cartesian(ylim = c(ymin, ymax)) +
ggtitle(paste0("ROI:", region))
plotobject
}plot_univar("superiortemporal_L", -2, 3)## Warning: Using alpha for a discrete variable is not advised.
plot_univar("superiortemporal_R", -2, 4)## Warning: Using alpha for a discrete variable is not advised.
plot_univar("supramarginal_L", -2, 4)## Warning: Using alpha for a discrete variable is not advised.
plot_univar("supramarginal_R", -2, 4)## Warning: Using alpha for a discrete variable is not advised.
plot_univar("V1_bilateral", -1, 11)## Warning: Using alpha for a discrete variable is not advised.
plot_univar("MT_bilateral", -1, 7)## Warning: Using alpha for a discrete variable is not advised.
plot_univar("antParietal_bilateral", -1, 7)## Warning: Using alpha for a discrete variable is not advised.
plot_univar("precentral_inferiorfrontal_R", -1, 7)## Warning: Using alpha for a discrete variable is not advised.
2.2 Check whether we should restrict to first 2 runs
Following the preregistration plan, we check whether we are able to take the second level cope values (averaging across 4 runs) as the main dependent measure, or whether we need to restrict the analysis to just the first 2 runs.
First we fit a model including an interaction between run number and event.
focal_data_checking <- focal_data %>%
filter(event != "fam",
top_n_voxels == "100",
domain != "both",
str_length(extracted_run_number) == 4) %>%
filter(ROI_name_final != "V1_bilateral",
ROI_name_final != "MT_bilateral")
focal_data_checking$fixation_position <- as.factor(focal_data_checking$fixation_position)
focal_data_checking$event <- relevel(focal_data_checking$event, ref = "unexp")
checkmodel <- lmer(
data = focal_data_checking,
formula = meanbeta ~ extracted_run_number * event + (1|subjectID) + (1|ROI_name_final)
)
checkmodel %>%
apa_print() %>%
apa_table()
**
| Term | \(\hat{\beta}\) | 95% CI | \(t\) | \(\mathit{df}\) | \(p\) |
|---|---|---|---|---|---|
| Intercept | 2.20 | [1.25, 3.16] | 4.51 | 6.96 | .003 |
| Extracted run number1 | 0.41 | [0.27, 0.55] | 5.83 | 3,004.27 | < .001 |
| Extracted run number2 | 0.15 | [0.02, 0.29] | 2.18 | 3,004.27 | .030 |
| Extracted run number3 | -0.15 | [-0.29, -0.01] | -2.06 | 3,006.55 | .039 |
| Event1 | 0.09 | [0.01, 0.17] | 2.27 | 3,003.97 | .023 |
| Extracted run number1 \(\times\) Event1 | 0.06 | [-0.08, 0.19] | 0.78 | 3,003.97 | .433 |
| Extracted run number2 \(\times\) Event1 | 0.07 | [-0.07, 0.21] | 1.00 | 3,003.97 | .318 |
| Extracted run number3 \(\times\) Event1 | -0.07 | [-0.21, 0.07] | -0.92 | 3,003.97 | .356 |
summary(checkmodel)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ extracted_run_number * event + (1 | subjectID) + (1 |
## ROI_name_final)
## Data: focal_data_checking
##
## REML criterion at convergence: 13783
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.973 -0.603 -0.053 0.550 5.440
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.19 1.09
## ROI_name_final (Intercept) 1.20 1.10
## Residual 5.14 2.27
## Number of obs: 3048, groups: subjectID, 32; ROI_name_final, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.2045 0.4890 6.9622 4.51 0.0028
## extracted_run_number1 0.4138 0.0710 3004.2738 5.83 6.1e-09
## extracted_run_number2 0.1544 0.0710 3004.2738 2.18 0.0297
## extracted_run_number3 -0.1484 0.0720 3006.5502 -2.06 0.0393
## event1 0.0934 0.0411 3003.9740 2.27 0.0230
## extracted_run_number1:event1 0.0556 0.0709 3003.9740 0.78 0.4332
## extracted_run_number2:event1 0.0709 0.0709 3003.9740 1.00 0.3180
## extracted_run_number3:event1 -0.0662 0.0717 3003.9740 -0.92 0.3557
##
## (Intercept) **
## extracted_run_number1 ***
## extracted_run_number2 *
## extracted_run_number3 *
## event1 *
## extracted_run_number1:event1
## extracted_run_number2:event1
## extracted_run_number3:event1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ext__1 ext__2 ext__3 event1 e__1:1 e__2:1
## extrctd_r_1 -0.001
## extrctd_r_2 -0.001 -0.329
## extrctd_r_3 0.002 -0.338 -0.338
## event1 0.000 0.000 0.000 0.000
## extrct__1:1 0.000 0.000 0.000 0.000 -0.005
## extrct__2:1 0.000 0.000 0.000 0.000 -0.005 -0.330
## extrct__3:1 0.000 0.000 0.000 0.000 0.014 -0.337 -0.337lsmeans(checkmodel, pairwise ~ event | extracted_run_number)$contrasts## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3048' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3048)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 3048' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3048)' or larger];
## but be warned that this may result in large computation time and memory use.## extracted_run_number = run1:
## contrast estimate SE df z.ratio p.value
## unexp - exp 0.298 0.164 Inf 1.821 0.0690
##
## extracted_run_number = run2:
## contrast estimate SE df z.ratio p.value
## unexp - exp 0.328 0.164 Inf 2.008 0.0450
##
## extracted_run_number = run3:
## contrast estimate SE df z.ratio p.value
## unexp - exp 0.054 0.166 Inf 0.327 0.7440
##
## extracted_run_number = run4:
## contrast estimate SE df z.ratio p.value
## unexp - exp 0.066 0.164 Inf 0.405 0.6850
##
## Degrees-of-freedom method: asymptoticmodel_habituation_results <- cbind(gen.m(checkmodel), gen.ci(checkmodel)[3:10,])## Computing profile confidence intervals ...model_habituation_results_byrun <- lsmeans(checkmodel, pairwise ~ event | extracted_run_number, pbkrtest.limit = 3048)$contrasts %>% as.data.frame()We then compare this model to a simpler model without this interaction.
checkmodel_simpler <- lmer(
data = focal_data_checking,
formula = meanbeta ~ event + extracted_run_number + (1|subjectID) + (1|ROI_name_final)
)
summary(checkmodel_simpler)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event + extracted_run_number + (1 | subjectID) + (1 |
## ROI_name_final)
## Data: focal_data_checking
##
## REML criterion at convergence: 13774
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.947 -0.602 -0.057 0.553 5.414
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.19 1.09
## ROI_name_final (Intercept) 1.20 1.10
## Residual 5.14 2.27
## Number of obs: 3048, groups: subjectID, 32; ROI_name_final, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.2045 0.4890 6.9620 4.51 0.0028 **
## event1 0.0939 0.0411 3006.9741 2.29 0.0222 *
## extracted_run_number1 0.4138 0.0710 3007.2741 5.83 6.1e-09 ***
## extracted_run_number2 0.1544 0.0710 3007.2741 2.18 0.0297 *
## extracted_run_number3 -0.1484 0.0720 3009.5523 -2.06 0.0393 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 ext__1 ext__2
## event1 0.000
## extrctd_r_1 -0.001 0.000
## extrctd_r_2 -0.001 0.000 -0.329
## extrctd_r_3 0.002 0.000 -0.338 -0.338# checkmodel_simpler %>%
# apa_print() %>%
# apa_table()
anova(checkmodel, checkmodel_simpler)## refitting model(s) with ML (instead of REML)knitr::kable(model_habituation_results)| B | SE | df | t | p | CI_95_lower | CI_95_upper | |
|---|---|---|---|---|---|---|---|
| (Intercept) | 2.204 | 0.489 | 6.96 | 4.508 | 0.003 | 2.208 | 2.323 |
| extracted_run_number1 | 0.414 | 0.071 | 3004.27 | 5.831 | 0.000 | 1.183 | 3.226 |
| extracted_run_number2 | 0.154 | 0.071 | 3004.27 | 2.175 | 0.030 | 0.275 | 0.553 |
| extracted_run_number3 | -0.148 | 0.072 | 3006.55 | -2.062 | 0.039 | 0.015 | 0.293 |
| event1 | 0.093 | 0.041 | 3003.97 | 2.274 | 0.023 | -0.289 | -0.007 |
| extracted_run_number1:event1 | 0.056 | 0.071 | 3003.97 | 0.784 | 0.433 | 0.013 | 0.174 |
| extracted_run_number2:event1 | 0.071 | 0.071 | 3003.97 | 0.999 | 0.318 | -0.083 | 0.195 |
| extracted_run_number3:event1 | -0.066 | 0.072 | 3003.97 | -0.924 | 0.356 | -0.068 | 0.210 |
knitr::kable(model_habituation_results_byrun)| contrast | extracted_run_number | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| unexp - exp | run1 | 0.298 | 0.164 | 3004 | 1.821 | 0.069 |
| unexp - exp | run2 | 0.328 | 0.164 | 3004 | 2.008 | 0.045 |
| unexp - exp | run3 | 0.054 | 0.166 | 3004 | 0.327 | 0.744 |
| unexp - exp | run4 | 0.066 | 0.164 | 3004 | 0.405 | 0.685 |
# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_perrun.Rds"))This is a grey area. On the one hand we do not observe a significant run x event interaction. On the other hand, the primary VOE effect (unexp > expected) is indeed variable across runs, and stronger in runs 1 and 2 than runs 3 and 4, when considering psychology and physics events). So on balance, we have the grounds for focusing the analysis on runs 1-2.
Here we fit a model including main effects of domain and event (domain-general and domain-specific regions) and also their interaction (domain-specific regions only), and organize and print the outputs.
2.3 Model Tables
data_for_modeling <- univariate_data %>%
filter(event != "fam",
# top_n_voxels == "100",
!str_detect(extracted_run_number, "fold"),
!str_detect(extracted_run_number, "all_runs"),
domain != "both")
ROIs <- levels(as.factor(data_for_modeling$ROI_name_final))
runs <- c("run1", "run2", "run3", "run4", "all_runs", "runs12")
top_n_voxels <- unique(data_for_modeling$top_n_voxels)
# top_n_voxels <- 100
effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)
rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "domain1", "event:domain")modelsummaries_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (topN in top_n_voxels) {
cur_top_voxels_N <- topN
for (run in runs) {
curRun = run
if (str_length(curRun) == 4) { # run1, run2, etc
ROIdata <- data_for_modeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN,
extracted_run_number == curRun)
} else if (curRun == "all_runs") {
ROIdata <- data_for_modeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN,
str_length(extracted_run_number) == 4)
} else if (curRun == "runs12") {
ROIdata <- data_for_modeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN) %>%
filter(extracted_run_number == "run1" | extracted_run_number == "run2")
}
# print("Current ROI: ", curROI)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
} else {
ROI_category <- "general"
}
# had to remove run random int (convergence))
model <- lmer(data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID))
model_just_maineffects <-
update(model, formula = ~ . - event:domain) # Without interaction term
model_just_domain <-
update(model_just_maineffects, formula = ~ . - event)
model_just_event <-
update(model_just_maineffects, formula = ~ . - domain)
BF_BIC_interaction <-
data.frame(BF = exp((
BIC(model_just_maineffects) - BIC(model)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_maineffects) - BIC(model)
) / 2))) # BICs to Bayes factor
BF_BIC_event <-
data.frame(BF = exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2)))
BF_BIC_domain <-
data.frame(BF = exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2)))
BF_initial <- rbind(BF_BIC_event,
BF_BIC_domain,
BF_BIC_interaction)
rownames(BF_initial) <- rownames_no_intercept
BFs <- rbind(BF_intercept_row,
BF_initial)
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <-
rbind(effect_sizes_intercept_row, effect_sizes_initial)
ROI_focal <- ROIdata$focal_region[1]
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:6, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d",
"BF",
"BF_interpretation"
)
modelsummaries_init <-
rbind(modelsummaries_init, modelsummary)
cat(
tab_model(
model,
show.std = TRUE,
show.se = TRUE,
show.stat = TRUE,
auto.label = TRUE,
title = paste0("Confirmatory Univariate ROI results in: ", curROI, " ", curRun, " top voxels = ", cur_top_voxels_N)
)$knitr,
" \n \n"
)
}
}
}modelsummaries <- modelsummaries_init %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
modelsummaries$effect <-
as.factor(modelsummaries$effect)
levels(modelsummaries$effect) <-
c("intercept", "domain", "event", "event:domain")
rownames(modelsummaries) <- NULL
# write_rds(modelsummaries, here("outputs/univariate_results/univar_results_alltopN_allregions.Rds"))modelsummaries <- read_rds(here("outputs/univariate_results/univar_results_alltopN_allregions.Rds"))DT::datatable(modelsummaries %>% filter(focal_region == 1, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect), options = list(scrollX = TRUE))Below we check, for each domain-specific region (left and right STS and SMG), whether a model including an interaction between domain and event fits better than a model without the interaction.
focal_data_100_runs12 <- focal_data %>%
filter(top_n_voxels == "100",
domain != "both",
event != "fam") %>%
filter(extracted_run_number == "run1" | extracted_run_number == "run2")LSMG_model <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
formula = meanbeta ~ event * domain + (1|subjectID)
)
plot(allEffects(LSMG_model))
# check_model(LSMG_model)
check_outliers(LSMG_model)## OK: No outliers detected.summary(LSMG_model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | subjectID)
## Data: focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L")
##
## REML criterion at convergence: 948
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4287 -0.5510 0.0345 0.5874 2.6185
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.02 2.00
## Residual 1.58 1.26
## Number of obs: 256, groups: subjectID, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.3899 0.3630 31.0000 6.58 2.4e-07 ***
## event1 0.0672 0.0787 221.0000 0.85 0.3942
## domain1 0.2580 0.0787 221.0000 3.28 0.0012 **
## event1:domain1 0.1877 0.0787 221.0000 2.39 0.0178 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1
## event1 0.000
## domain1 0.000 0.000
## event1:dmn1 0.000 0.000 0.000LSMG_results_bydomain_top100 <-
lsmeans(LSMG_model, pairwise ~ event |
domain)$contrasts %>% as.data.frame()
knitr::kable(LSMG_results_bydomain_top100)| contrast | domain | estimate | SE | df | t.ratio | p.value |
|---|---|---|---|---|---|---|
| unexp - exp | physics | 0.510 | 0.222 | 221 | 2.29 | 0.023 |
| unexp - exp | psychology | -0.241 | 0.222 | 221 | -1.08 | 0.280 |
# write_rds(LSMG_results_bydomain_top100, path = here("outputs/univariate_results/LSMG_results_bydomain_top100.Rds"))
LSMG_model_simpler <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
formula = meanbeta ~ event + domain + (1|subjectID)
)
# check_model(LSMG_model_simpler)
check_outliers(LSMG_model_simpler)## OK: No outliers detected.LSMG_model_comparison <- anova(LSMG_model_simpler, LSMG_model) %>%
as.data.frame()## refitting model(s) with ML (instead of REML)# write_rds(LSMG_model_comparison, path = here("outputs/univariate_results/LSMG_model_comparison_top100.Rds"))RSMG_model <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_R"),
formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
)
summary(RSMG_model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |
## subjectID)
## Data: focal_data_100_runs12 %>% filter(ROI_name_final == "supramarginal_R")
##
## REML criterion at convergence: 996
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.437 -0.585 0.036 0.515 3.808
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 6.484 2.546
## extracted_run_number (Intercept) 0.098 0.313
## Residual 1.829 1.353
## Number of obs: 256, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.5420 0.5087 14.6472 5.00 0.00017 ***
## event1 0.1613 0.0845 219.9999 1.91 0.05767 .
## domain1 0.5037 0.0845 219.9999 5.96 1e-08 ***
## event1:domain1 0.0566 0.0845 219.9999 0.67 0.50410
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1
## event1 0.000
## domain1 0.000 0.000
## event1:dmn1 0.000 0.000 0.000#
# RSMG_model %>%
# apa_print() %>%
# apa_table()
# check_model(RSMG_model)
check_outliers(RSMG_model)## OK: No outliers detected.lsmeans(RSMG_model, pairwise ~ event | domain)$contrasts## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.436 0.239 220 1.822 0.0700
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.209 0.239 220 0.876 0.3820
##
## Degrees-of-freedom method: kenward-rogerRSTS_model <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_R"),
formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
)
summary(RSTS_model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |
## subjectID)
## Data:
## focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_R")
##
## REML criterion at convergence: 960
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4296 -0.5228 0.0122 0.5633 3.0469
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 2.6191 1.618
## extracted_run_number (Intercept) 0.0135 0.116
## Residual 1.7596 1.326
## Number of obs: 256, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.0199 0.3090 20.7114 6.54 1.9e-06 ***
## event1 0.1234 0.0829 220.0000 1.49 0.1381
## domain1 -0.2421 0.0829 220.0000 -2.92 0.0039 **
## event1:domain1 0.0539 0.0829 220.0000 0.65 0.5160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1
## event1 0.000
## domain1 0.000 0.000
## event1:dmn1 0.000 0.000 0.000# RSTS_model %>%
# apa_print() %>%
# apa_table()
# check_model(RSTS_model)
check_outliers(RSTS_model)## OK: No outliers detected.lsmeans(RSTS_model, pairwise ~ event | domain)$contrasts## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.355 0.234 220 1.512 0.1320
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.139 0.234 220 0.592 0.5540
##
## Degrees-of-freedom method: kenward-rogerLSTS_model <- lmer(
data = focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_L"),
formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
)
summary(LSTS_model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + (1 | extracted_run_number) + (1 |
## subjectID)
## Data:
## focal_data_100_runs12 %>% filter(ROI_name_final == "superiortemporal_L")
##
## REML criterion at convergence: 1041
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.539 -0.593 0.035 0.503 3.454
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.337 2.083
## extracted_run_number (Intercept) 0.142 0.376
## Residual 2.362 1.537
## Number of obs: 256, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.6181 0.4644 6.5575 1.33 0.228
## event1 0.0624 0.0960 220.0000 0.65 0.517
## domain1 -0.2998 0.0960 220.0000 -3.12 0.002 **
## event1:domain1 -0.0658 0.0960 220.0000 -0.69 0.494
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1
## event1 0.000
## domain1 0.000 0.000
## event1:dmn1 0.000 0.000 0.000# LSTS_model %>%
# apa_print() %>%
# apa_table()
# check_model(LSTS_model)
check_outliers(LSTS_model)## Warning: 1 outliers detected (cases 52).outliers_list <- check_outliers(LSTS_model)
insight::get_data(LSTS_model)[outliers_list, ]LSTS_model_nooutliers <- lmer(
data = insight::get_data(LSTS_model) %>%
filter(meanbeta != -11.6),
formula = meanbeta ~ event * domain + (1|extracted_run_number) + (1|subjectID)
)
# summary(LSTS_model_nooutliers)
LSTS_model_nooutliers %>%
apa_print() %>%
apa_table()
**
| Term | \(\hat{\beta}\) | 95% CI | \(t\) | \(\mathit{df}\) | \(p\) |
|---|---|---|---|---|---|
| Intercept | 0.62 | [-0.29, 1.53] | 1.33 | 6.56 | .228 |
| Event1 | 0.06 | [-0.13, 0.25] | 0.65 | 220.00 | .517 |
| Domain1 | -0.30 | [-0.49, -0.11] | -3.12 | 220.00 | .002 |
| Event1 \(\times\) Domain1 | -0.07 | [-0.25, 0.12] | -0.69 | 220.00 | .494 |
lsmeans(LSTS_model_nooutliers, pairwise ~ event | domain)$contrasts## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp -0.0069 0.272 220 -0.026 0.9800
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.2564 0.272 220 0.944 0.3460
##
## Degrees-of-freedom method: kenward-roger3 PART 1b: Sensitivity to voxel N
focal_modelsummaries <- modelsummaries %>%
filter(focal_region == 1)
focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "runs12", "all_runs"))
ggplot(focal_modelsummaries %>%
filter(str_length(extracted_run_number) == 4,
effect != "intercept"), aes(extracted_run_number, d, colour = top_n_voxels)) +
geom_point() +
geom_line(aes(group = top_n_voxels)) +
ylab("Cohen's D") +
facet_wrap(~ effect + factor(ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_bilateral", "MT_bilateral")), nrow = 3) +
# coord_cartesian(ylim = c(-.75, .75)) +
geom_hline(yintercept = 0) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
4 PART 2: All regions (exploratory)
Then we repeat the same modeling procedure on these regions. But this time, we fit an interaction between event and domain in all regions, because some MD regions appear to have a preference for physical events, and it is an open question, if this is true, whether they too could encode physical prediction error.
DT::datatable(modelsummaries %>% filter(focal_region == 1, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect))This set of regions includes left and right V1, MT, insula, and IFG (rather than bilateral). The bilateral versions of these ROIs are not included in this portion of the analysis.
DT::datatable(modelsummaries %>% filter(focal_region == 0, top_n_voxels == 100, extracted_run_number == "runs12") %>% arrange(effect), options = list(scrollX = TRUE))5 PART 3: Action-env scenarios
5.1 Check for run X event effects
focal_data_both_100 <- focal_data %>%
filter(event != "fam",
top_n_voxels == "100",
domain == "both")options(contrasts = c("contr.sum", "contr.poly"))
checkdata_both <- focal_data_both_100 %>%
filter(ROI_name_final != "V1_bilateral",
ROI_name_final != "MT_bilateral",
str_length(extracted_run_number) == 4)
checkdata_both$fixation_position <- as.factor(checkdata_both$fixation_position)
checkdata_both$event <- relevel(checkdata_both$event, ref = "unexp")
checkmodel_bothdomains <- lmer(
data = checkdata_both,
formula = meanbeta ~ event * extracted_run_number + (1|subjectID) + (1|ROI_name_final)
)
summary(checkmodel_bothdomains)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * extracted_run_number + (1 | subjectID) + (1 |
## ROI_name_final)
## Data: checkdata_both
##
## REML criterion at convergence: 7354
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.428 -0.547 -0.028 0.538 5.262
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.11 1.05
## ROI_name_final (Intercept) 1.34 1.16
## Residual 6.83 2.61
## Number of obs: 1524, groups: subjectID, 32; ROI_name_final, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.0973 0.5121 6.6088 4.10 0.0052
## event1 0.2175 0.0669 1479.8794 3.25 0.0012
## extracted_run_number1 0.5830 0.1157 1480.2519 5.04 5.2e-07
## extracted_run_number2 0.1838 0.1157 1480.2519 1.59 0.1123
## extracted_run_number3 -0.4030 0.1172 1483.0707 -3.44 0.0006
## event1:extracted_run_number1 -0.1024 0.1156 1479.8794 -0.89 0.3760
## event1:extracted_run_number2 -0.0745 0.1156 1479.8794 -0.64 0.5194
## event1:extracted_run_number3 0.2402 0.1169 1479.8794 2.06 0.0400
##
## (Intercept) **
## event1 **
## extracted_run_number1 ***
## extracted_run_number2
## extracted_run_number3 ***
## event1:extracted_run_number1
## event1:extracted_run_number2
## event1:extracted_run_number3 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 ext__1 ext__2 ext__3 e1:__1 e1:__2
## event1 0.000
## extrctd_r_1 -0.001 0.000
## extrctd_r_2 -0.001 0.000 -0.329
## extrctd_r_3 0.002 0.000 -0.338 -0.338
## evnt1:xt__1 0.000 -0.005 0.000 0.000 0.000
## evnt1:xt__2 0.000 -0.005 0.000 0.000 0.000 -0.330
## evnt1:xt__3 0.000 0.014 0.000 0.000 0.000 -0.337 -0.337# checkmodel %>%
# apa_print() %>%
# apa_table()
plot(allEffects(checkmodel_bothdomains))
lsmeans(checkmodel_bothdomains, pairwise ~ event | extracted_run_number)$contrasts## extracted_run_number = run1:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.230 0.267 1480 0.860 0.3880
##
## extracted_run_number = run2:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.286 0.267 1480 1.070 0.2840
##
## extracted_run_number = run3:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.915 0.271 1480 3.380 0.0010
##
## extracted_run_number = run4:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.308 0.267 1480 1.160 0.2480
##
## Degrees-of-freedom method: kenward-rogermodel_habituation_results <- cbind(gen.m(checkmodel_bothdomains), gen.ci(checkmodel_bothdomains)[3:10,])## Computing profile confidence intervals ...model_habituation_results_byrun <- lsmeans(checkmodel_bothdomains, pairwise ~ event | extracted_run_number, pbkrtest.limit = 3048)$contrasts %>% as.data.frame()# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_psych-env_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_psych-env_perrun.Rds"))univariate_data_both <- univariate_data %>%
filter(event != "fam",
# top_n_voxels == "100",
domain == "both",
str_length(extracted_run_number) == 4)
ROIs <- levels(as.factor(univariate_data_both$ROI_name_final))
runs <- c("run1", "run2", "run3", "run4", "all_runs", "runs12")
# top_n_voxels <- 100
modelsummaries_both_init <- data.frame()5.2 Model tables
for (ROI in ROIs) {
curROI <- ROI
for (topN in top_n_voxels) {
cur_top_voxels_N <- topN
for (run in runs) {
curRun <- run
if (curRun != "all_runs" & curRun != "runs12") {
ROIdata <- univariate_data_both %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN,
extracted_run_number == curRun)
} else if (curRun == "all_runs") {
ROIdata <- univariate_data_both %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN)
} else if (curRun == "runs12") {
ROIdata <- univariate_data_both %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN) %>%
filter(extracted_run_number == "run1" |
extracted_run_number == "run2")
}
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
}
else {
ROI_category <- "general"
}
model <- lmer(data = ROIdata,
formula = meanbeta ~ event + (1 | subjectID))
focal_region <- ROIdata$focal_region[1]
model_null <-
update(model, formula = ~ . - event) # Without interaction term
BF_BIC_event <-
data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_null) - BIC(model)
) / 2))) # BICs to Bayes factor
# BF_initial <- rbind(BF_BIC_event)
rownames(BF_BIC_event) <- c("event1")
BFs <- rbind(BF_intercept_row,
BF_BIC_event)
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <-
rbind(effect_sizes_intercept_row, effect_sizes_initial)
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
focal_region,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:4,],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d",
"BF",
"BF_interpretation"
)
modelsummaries_both_init <-
rbind(modelsummaries_both_init, modelsummary)
cat(
tab_model(
model,
show.std = TRUE,
show.se = TRUE,
show.stat = TRUE,
auto.label = TRUE,
title = paste0(
"Exploratory Univariate ROI results for events from both domains, run = ",
curRun,
curROI
)
)$knitr,
" \n \n"
)
}
}
}modelsummaries_both <- modelsummaries_both_init %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
modelsummaries_both$effect <-
as.factor(modelsummaries_both$effect)
levels(modelsummaries_both$effect) <-
c("intercept", "event")
rownames(modelsummaries_both) <- NULL
# write_rds(modelsummaries_both, here("outputs/univariate_analysis_outputs/univar_results_psych-env_allregions_alltopN.Rds"))modelsummaries_both <- read_rds( here("outputs/univariate_results/univar_results_psych-env_allregions_alltopN.Rds"))5.3 Visualize
Effects by run and ROI size
modelsummaries_both$extracted_run_number <- factor(modelsummaries_both$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "runs12", "all_runs"))
modelsummaries_both_focal <- modelsummaries_both %>% filter(focal_region == 1, effect == "event")
ggplot(modelsummaries_both_focal %>% filter(str_length(extracted_run_number) == 4), aes(extracted_run_number, d, colour = top_n_voxels)) +
geom_point() +
geom_line(aes(group = top_n_voxels)) +
# geom_errorbar(aes(ymin = B - SE, ymax = B + SE, width = 0.1)) +
# ylab("Bayes Factor (log)") +
ylab("Cohen's D") +
facet_wrap(~ effect + ROI, nrow = 1) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
5.3.1 Snsitivity to top voxel N
focal_modelsummaries <- modelsummaries_both %>%
filter(focal_region == 1,
str_length(extracted_run_number) == 4)
focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4"))
ggplot(
focal_modelsummaries %>%
filter(str_length(extracted_run_number) == 4,
effect == "event"),
aes(extracted_run_number, d, colour = top_n_voxels)
) +
geom_point() +
geom_line(aes(group = top_n_voxels)) +
ylab("Cohen's D") +
facet_wrap( ~ effect + factor(
ROI,
levels = c(
"supramarginal_L",
"supramarginal_R",
"superiortemporal_L",
"superiortemporal_R",
"antParietal_bilateral",
"precentral_inferiorfrontal_R",
"V1_bilateral",
"MT_bilateral"
)
), nrow = 1) +
# coord_cartesian(ylim = c(-.75, .75)) +
geom_hline(yintercept = 0) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) 
5.4 All regions, all runs
DT::datatable(modelsummaries_both %>% filter(extracted_run_number == "all_runs") %>% arrange(effect, focal_region), options = list(scrollX = TRUE))5.5 All regions, runs 1-2
DT::datatable(modelsummaries_both %>% filter(extracted_run_number == "runs12") %>% arrange(effect, focal_region), options = list(scrollX = TRUE))6 PART 4: Using VOE runs 3-4 to localize
Here we repeat the same analysis as confirmatory analyses above, except that we use VOE runs 3-4 to localize STS and SMG, and we search for voxels in the domain general ROI parcels that respond more to unexp than exp, to ask whether there are any voxels that show the VOE response in those regions, without considering the localizer data.
univariate_data_VOEextract <-
readRDS(here("outputs", "univariate_data", "study2_univariate_data_VOEextract.Rds"))
univariate_summary_domain_VOEextract <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_VOEextract.Rds"))
univariate_summary_domain_runs12_VOEextract <-
readRDS(here("outputs", "descriptive_summaries", "study2_univariate_summary_domain_runs12_VOEextract.Rds"))focal_data_VOEextract <- univariate_data_VOEextract %>%
filter(focal_region == 1)
focal_summary_domain_VOEextract <- univariate_summary_domain_VOEextract %>%
filter(focal_region == 1)
focal_summary_domain_runs12_VOEextract <- univariate_summary_domain_runs12_VOEextract %>%
filter(focal_region == 1)
domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")
regions <- levels(as.factor(focal_data_VOEextract$ROI_name_final))6.1 Plots
First, here are plots of betas, run by run, to physics and psychology events (familiarization, expected, unexpected).
plot_univar_runbyrun("superiortemporal_L")## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("superiortemporal_R")## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("supramarginal_L")## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("supramarginal_R")## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("antParietal_bilateral")## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("precentral_inferiorfrontal_R")## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("V1_bilateral")## Warning: Using alpha for a discrete variable is not advised.
plot_univar_runbyrun("MT_bilateral")## Warning: Using alpha for a discrete variable is not advised.
DS_runs12_plots <- plot_univar_runs12("superiortemporal", 0, 4) |
plot_univar_runs12("supramarginal", 0, 4)
DG_runs12_plots <- plot_univar_runs12("precentral_inferiorfrontal_R", 0, 4) | plot_univar_runs12("antParietal_bilateral", 0, 4) | plot_univar_runs12("V1_bilateral", 0, 11) | plot_univar_runs12("MT_bilateral", 0, 11)
DS_runs12_plots / DG_runs12_plots## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
focal_data_100_runs12_VOEextract <- focal_data_VOEextract %>%
filter(
event != "fam",
top_n_voxels == "100",
domain != "both",
(extracted_run_number == "run1" | extracted_run_number == "run2"))
ROIs <- levels(as.factor(focal_data_100_runs12_VOEextract$ROI_name_final))
focal_data_100_runs12_VOEextract$event <- relevel(focal_data_100_runs12_VOEextract$event, ref = "exp")
focal_modelsummaries_100_VOEextract_init <- data.frame()6.2 Model Tables
6.2.1 Psych-action and physics events
for (ROI in ROIs) {
curROI <- ROI
ROIdata <-
focal_data_100_runs12_VOEextract %>% filter(ROI_name_final == curROI)
# print("Current ROI: ", curROI)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
model <- lmer(
data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID)) # had to remove run random int (convergence))
} else {
ROI_category <- "general"
model <- lmer(data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID))
}
model_just_maineffects <- update(model, formula = ~ . -event:domain) # Without interaction term
model_just_domain <- update(model_just_maineffects, formula = ~ . -event)
model_just_event <- update(model_just_maineffects, formula = ~ . -domain)
BF_BIC_interaction <- data.frame(BF = exp((BIC(model_just_maineffects) - BIC(model))/2), BF_interpretation = interpret_bf(exp((BIC(model_just_maineffects) - BIC(model))/2))) # BICs to Bayes factor
BF_BIC_event <- data.frame(BF = exp((BIC(model_just_domain) - BIC(model_just_maineffects))/2),
BF_interpretation = interpret_bf(exp((BIC(model_just_domain) - BIC(model_just_maineffects))/2)))
BF_BIC_domain <- data.frame(BF = exp((BIC(model_just_event) - BIC(model_just_maineffects))/2),
BF_interpretation = interpret_bf(exp((BIC(model_just_event) - BIC(model_just_maineffects))/2)))
BF_initial <- rbind(BF_BIC_event,
BF_BIC_domain,
BF_BIC_interaction)
rownames(BF_initial) <- rownames_no_intercept
BFs <- rbind(BF_intercept_row,
BF_initial)
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <- rbind(effect_sizes_intercept_row, effect_sizes_initial)
ROI_focal <- ROIdata$focal_region[1]
modelsummary <-
cbind(
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:6, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d",
"BF",
"BF_interpretation"
)
focal_modelsummaries_100_VOEextract_init <-
rbind(focal_modelsummaries_100_VOEextract_init, modelsummary)
# model_parameters(model) %>% print_md()
cat(
tab_model(
model,
show.std = TRUE,
show.se = TRUE,
show.stat = TRUE,
auto.label = TRUE,
title = paste0("Explratory Univariate ROI results (localized using VOE runs 3-4) in: ", curROI)
)$knitr,
" \n \n"
)
# focal_modeltables_100 <- append(focal_modeltables_100, cur_tab)
}focal_modelsummaries_100_VOEextract <- focal_modelsummaries_100_VOEextract_init %>%
mutate_at(c(5:13), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
focal_modelsummaries_100_VOEextract$effect <-
as.factor(focal_modelsummaries_100_VOEextract$effect)
levels(focal_modelsummaries_100_VOEextract$effect) <-
c("intercept", "domain", "event", "event:domain")
rownames(focal_modelsummaries_100_VOEextract) <- NULL
# write_rds(focal_modelsummaries_100_VOEextract, here("outputs/univariate_results/univar_results_top100_VOE-extract_runs12.Rds"))focal_modelsummaries_100_VOEextract <- read_rds(here("outputs/univariate_results/univar_results_top100_VOE-extract_runs12.Rds"))
DT::datatable(dplyr::arrange(focal_modelsummaries_100_VOEextract, effect), options = list(scrollX = TRUE))6.2.2 Psych-env events
focal_data_100_runs1234_VOEextract_bothdomains <-
focal_data_VOEextract %>%
filter(
event != "fam",
top_n_voxels == "100",
domain == "both",
str_length(extracted_run_number) == 4
)
ROIs <-
levels(as.factor(
focal_data_100_runs1234_VOEextract_bothdomains$ROI_name_final
))
focal_data_100_runs1234_VOEextract_bothdomains$event <-
relevel(focal_data_100_runs1234_VOEextract_bothdomains$event, ref = "unexp")
focal_modelsummaries_100_VOEextract_bothdomains_init <- data.frame()for (ROI in ROIs) {
curROI <- ROI
cur_top_voxels_N <- 100
curRun <- "runs12"
ROIdata <- focal_data_100_runs1234_VOEextract_bothdomains %>%
filter(ROI_name_final == curROI) %>%
filter(extracted_run_number == "run1" |
extracted_run_number == "run2")
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
}
else {
ROI_category <- "general"
}
model <- lmer(data = ROIdata,
formula = meanbeta ~ event + (1 | subjectID))
focal_region <- ROIdata$focal_region[1]
model_null <-
update(model, formula = ~ . - event) # Without interaction term
BF_BIC_event <-
data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_null) - BIC(model)
) / 2))) # BICs to Bayes factor
# BF_initial <- rbind(BF_BIC_event)
rownames(BF_BIC_event) <- c("event1")
BFs <- rbind(BF_intercept_row,
BF_BIC_event)
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <-
rbind(effect_sizes_intercept_row, effect_sizes_initial)
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
focal_region,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:4, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d",
"BF",
"BF_interpretation"
)
focal_modelsummaries_100_VOEextract_bothdomains_init <-
rbind(focal_modelsummaries_100_VOEextract_bothdomains_init,
modelsummary)
cat(
tab_model(
model,
show.std = TRUE,
show.se = TRUE,
show.stat = TRUE,
auto.label = TRUE,
title = paste0(
"Exploratory Univariate ROI results for events from both domains, run = ",
curRun,
curROI
)
)$knitr,
" \n \n"
)
}focal_modelsummaries_100_VOEextract_bothdomains <- focal_modelsummaries_100_VOEextract_bothdomains_init %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
focal_modelsummaries_100_VOEextract_bothdomains$effect <-
as.factor(focal_modelsummaries_100_VOEextract_bothdomains$effect)
levels(focal_modelsummaries_100_VOEextract_bothdomains$effect) <-
c("intercept", "event")
rownames(focal_modelsummaries_100_VOEextract_bothdomains) <- NULL
# write_rds(focal_modelsummaries_100_VOEextract_bothdomains, here("outputs/univariate_results/univar_results_psych-env_top100_VOE-extract_runs12.Rds"))focal_modelsummaries_100_VOEextract_bothdomains <- read_rds( here("outputs/univariate_results/univar_results_psych-env_top100_VOE-extract_runs12.Rds"))DT::datatable(focal_modelsummaries_100_VOEextract_bothdomains, options = list(scrollX = TRUE))7 PART 5: Single betas
single_betas0 <- readRDS(here("outputs/univariate_data/study2_univariate_data_singlebetas.Rds")) %>%
select(-ROI_category)
regioninfo_category <- regioninfo %>%
select(ROI_name_final, ROI_category)
single_betas <- left_join(single_betas0, regioninfo_category)## Joining, by = "ROI_name_final"single_betas$ROI_name_final <- as.factor(single_betas$ROI_name_final)
single_betas$event <- as.factor(single_betas$event)
single_betas$event <- relevel(single_betas$event, ref = "fam")
single_betas$task <- factor(single_betas$task, levels = c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint"))single_betas_summary_domain <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_domain_singlebetas.Rds"))
single_betas_summary_task <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_task_singlebetas.Rds"))
single_betas_summary_scenario <- readRDS(here("outputs/descriptive_summaries/study2_univariate_summary_scenario_singlebetas.Rds"))
single_betas_summary_scenario$event <- relevel(single_betas_summary_scenario$event, ref = "fam")single_betas_focal <- single_betas %>%
filter(focal_region == 1)7.1 Plots
# solidity - 25a88b
# permanence - 19626a
# goal - f24b00
# efficiency - deae1e
# agent-solidity - 6600cd
# missing-barrier - fba0d0
ggplot(
single_betas_focal %>%
filter(
str_length(extracted_run_number) == 4,
event != "fam"
),
aes(task, meanbeta, alpha = event, fill = task)
) +
geom_bar(
position = "dodge",
stat = "summary",
fun.data = "mean_se",
colour = "black"
) +
geom_errorbar(
position = "dodge",
stat = "summary",
fun.data = "mean_se",
colour = "black"
) +
scale_fill_manual(values = c(
"#25a88b",
"#19626a",
"#f24b00",
"#deae1e",
"#6600cd",
"#fba0d0"
)) + # ,
# geom_point() +
theme(axis.text.x = element_blank()) +
# scale_alpha_continuous(range = c(0.5, 1)) +
facet_wrap( ~ ROI_name_final, scales = "free", nrow = 2) +
ggtitle("Runs 1-4, Focal regions")## Warning: Using alpha for a discrete variable is not advised.
ggplot(
single_betas_focal %>%
filter(ROI_name_final == "supramarginal_L"),
aes(
vis_statistics_video_group,
meanbeta,
alpha = event,
fill = task
)
) +
geom_bar(
position = "dodge",
stat = "summary",
fun.data = "mean_se",
colour = "black"
) +
geom_errorbar(
position = "dodge",
stat = "summary",
fun.data = "mean_se",
colour = "black"
) +
scale_fill_manual(values = c(
"#25a88b",
"#19626a",
"#f24b00",
"#deae1e",
"#6600cd",
"#fba0d0"
)) +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
ggtitle("Region: supramarginal_L")## Warning: Using alpha for a discrete variable is not advised.
Plotting each event type in each task against run number, per ROI.
ggplot(single_betas_focal, aes(as.integer(extracted_run_number), meanbeta, colour = event)) +
geom_smooth(method = "glm") +
facet_wrap( ~ ROI_name_final + task, nrow = 3, scales = "free") +
scale_colour_brewer(palette = "Set2")## `geom_smooth()` using formula = 'y ~ x'
Visualizing visual features against responses.
ggplot(single_betas_focal,
aes(normalized_iqr_contrast, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_luminance, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_motion, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_hsf, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_lsf, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_rect, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
ggplot(single_betas_focal,
aes(normalized_mean_curv, meanbeta)) +
geom_point() +
geom_smooth(method = "lm") +
facet_wrap( ~ ROI_name_final, nrow = 2, scales = "free")## `geom_smooth()` using formula = 'y ~ x'
7.2 Modeling
single_betas_focal$event <- relevel(single_betas_focal$event, ref = "unexp")single_betas_focal_psychphys_runs12 <- single_betas_focal %>%
filter(domain != "psychology-environment",
event != "fam") %>%
filter(extracted_run_number == "run1" | extracted_run_number == "run2")
ROIs <- unique(single_betas_focal_psychphys_runs12$ROI_name_final) %>% as.vector()
single_betas_focal_psychphys_runs12_summaries <- data.frame()
rownames_no_intercept <- c("event1", "domain1", "event:domain")7.3 Model Tables
for (ROI in ROIs) {
curROI <- ROI
ROIdata <-
single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == curROI)
# print("Current ROI: ", curROI)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
model <- lmer(
data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID) + (1|vis_statistics_video_group)
) # had to remove run random int (convergence))
} else {
ROI_category <- "general"
model <- lmer(
data = ROIdata,
formula = meanbeta ~ event * domain + (1|subjectID) + (1|vis_statistics_video_group)
)
}
model_just_maineffects <-
update(model, formula = ~ . - event:domain) # Without interaction term
model_just_domain <-
update(model_just_maineffects, formula = ~ . - event)
model_just_event <-
update(model_just_maineffects, formula = ~ . - domain)
BF_BIC_interaction <-
data.frame(BF = exp((
BIC(model_just_maineffects) - BIC(model)
) / 2), BF_interpretation = interpret_bf(exp((
BIC(model_just_maineffects) - BIC(model)
) / 2))) # BICs to Bayes factor
BF_BIC_event <-
data.frame(BF = exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2)))
BF_BIC_domain <-
data.frame(BF = exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2)))
BF_initial <- rbind(BF_BIC_event,
BF_BIC_domain,
BF_BIC_interaction)
rownames(BF_initial) <- rownames_no_intercept
BFs <- rbind(BF_intercept_row,
BF_initial)
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <-
rbind(effect_sizes_intercept_row, effect_sizes_initial)
ROI_focal <- ROIdata$focal_region[1]
cur_top_voxels_N <- 100
curRun <- "runs12"
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:6,],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d",
"BF",
"BF_interpretation"
)
single_betas_focal_psychphys_runs12_summaries <-
rbind(single_betas_focal_psychphys_runs12_summaries,
modelsummary)
# model_parameters(model) %>% print_md()
cat(
tab_model(
model,
show.std = TRUE,
show.se = TRUE,
show.stat = TRUE,
auto.label = TRUE,
title = paste0("Single beta Univariate ROI results (runs 1-2) in: ", curROI)
)$knitr,
" \n \n"
)
# focal_modeltables_100 <- append(focal_modeltables_100, cur_tab)
}single_betas_focal_psychphys_runs12_summaries <- single_betas_focal_psychphys_runs12_summaries %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
single_betas_focal_psychphys_runs12_summaries$effect <-
as.factor(single_betas_focal_psychphys_runs12_summaries$effect)
levels(single_betas_focal_psychphys_runs12_summaries$effect) <-
c("intercept", "domain", "event", "event:domain")
rownames(single_betas_focal_psychphys_runs12_summaries) <- NULL
# write_rds(single_betas_focal_psychphys_runs12_summaries, here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys.Rds"))single_betas_focal_psychphys_runs12_summaries <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys.Rds"))
DT::datatable(dplyr::arrange(single_betas_focal_psychphys_runs12_summaries, effect), options = list(scrollX = TRUE))7.4 PART 5a: Visual statistics (psych and physics)
single_betas_focal_visualstats_summary_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
ROIdata <-
single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == curROI)
# print("Current ROI: ", curROI)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "specific"
} else {
ROI_category <- "general"
}
model <- lmer(data = ROIdata, formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <-
rbind(effect_sizes_intercept_row, effect_sizes_initial)
ROI_focal <- ROIdata$focal_region[1]
cur_top_voxels_N <- 100
curRun <- "runs12"
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:12,],
effect_sizes
# BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d"
)
single_betas_focal_visualstats_summary_init <-
rbind(single_betas_focal_visualstats_summary_init,
modelsummary)
# model_parameters(model) %>% print_md()
# cat(
# tab_model(
# model,
# show.std = TRUE,
# show.se = TRUE,
# show.stat = TRUE,
# auto.label = TRUE,
# title = paste0("Exploratory Univariate ROI results in: ", curROI)
# )$knitr,
# " \n \n"
# )
}single_betas_focal_visualstats_summary <- single_betas_focal_visualstats_summary_init %>%
mutate_at(c(7:14), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
single_betas_focal_visualstats_summary$effect <-
as.factor(single_betas_focal_visualstats_summary$effect)
levels(single_betas_focal_visualstats_summary$effect)[1:4] <-
c("intercept", "domain", "event", "event:domain")
rownames(single_betas_focal_visualstats_summary) <- NULL
# write_rds(single_betas_focal_visualstats_summary, here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys_visstats.Rds"))single_betas_focal_visualstats_summary <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_psych_phys_visstats.Rds"))
DT::datatable(dplyr::arrange(single_betas_focal_visualstats_summary, effect), options = list(scrollX = TRUE))V1_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "V1_bilateral"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(V1_model_singlebetas)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "V1_bilateral")
##
## REML criterion at convergence: 3916
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.065 -0.581 -0.022 0.648 4.223
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 12.5313 3.540
## extracted_run_number (Intercept) 0.0802 0.283
## Residual 8.1251 2.850
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.0972 0.6683 25.6178 13.61 3.1e-13 ***
## event1 0.0655 0.1033 726.0307 0.63 0.5259
## domain1 0.4682 0.1767 726.0014 2.65 0.0082 **
## normalized_iqr_contrast -0.0388 0.1052 726.0004 -0.37 0.7119
## normalized_mean_luminance -1.2349 0.2240 726.0008 -5.51 4.9e-08 ***
## normalized_mean_motion 0.7183 0.1649 726.0038 4.36 1.5e-05 ***
## normalized_mean_hsf 0.3798 0.2439 726.0013 1.56 0.1198
## normalized_mean_rect -0.3466 0.2017 726.0002 -1.72 0.0861 .
## normalized_mean_curv 0.1841 0.1807 726.0021 1.02 0.3086
## event1:domain1 0.0613 0.1030 726.0002 0.60 0.5520
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.003
## domain1 -0.003 -0.055
## nrmlzd_qr_c -0.062 -0.036 0.084
## nrmlzd_mn_l 0.048 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.040 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.005 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.016 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.042 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.002 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023check_collinearity(V1_model_singlebetas)MT_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "MT_bilateral"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(MT_model_singlebetas)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "MT_bilateral")
##
## REML criterion at convergence: 3098
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.229 -0.591 0.009 0.605 3.287
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 6.351 2.520
## extracted_run_number (Intercept) 0.016 0.126
## Residual 2.716 1.648
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.42973 0.45984 30.85320 11.81 5.6e-13 ***
## event1 0.04583 0.05970 726.04405 0.77 0.44295
## domain1 0.18927 0.10213 726.00118 1.85 0.06427 .
## normalized_iqr_contrast -0.53523 0.06080 725.99966 -8.80 < 2e-16 ***
## normalized_mean_luminance -0.25298 0.12949 726.00030 -1.95 0.05112 .
## normalized_mean_motion 0.35157 0.09535 726.00465 3.69 0.00024 ***
## normalized_mean_hsf -1.38930 0.14100 726.00108 -9.85 < 2e-16 ***
## normalized_mean_rect 0.23006 0.11660 725.99933 1.97 0.04887 *
## normalized_mean_curv -0.59928 0.10449 726.00211 -5.74 1.4e-08 ***
## event1:domain1 0.00941 0.05954 725.99932 0.16 0.87444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.002
## domain1 -0.003 -0.055
## nrmlzd_qr_c -0.052 -0.036 0.084
## nrmlzd_mn_l 0.040 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.034 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.004 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.013 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.035 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.002 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023check_collinearity(MT_model_singlebetas)APC_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "antParietal_bilateral"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(APC_model_singlebetas)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "antParietal_bilateral")
##
## REML criterion at convergence: 3635
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.046 -0.535 0.015 0.520 4.165
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.3797 2.093
## extracted_run_number (Intercept) 0.0278 0.167
## Residual 5.7731 2.403
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.7034 0.4017 23.0478 9.22 3.4e-09 ***
## event1 0.2533 0.0870 726.0513 2.91 0.00372 **
## domain1 0.1349 0.1489 726.0023 0.91 0.36517
## normalized_iqr_contrast -0.2247 0.0887 726.0005 -2.54 0.01145 *
## normalized_mean_luminance -0.3695 0.1888 726.0013 -1.96 0.05073 .
## normalized_mean_motion 0.2846 0.1390 726.0063 2.05 0.04095 *
## normalized_mean_hsf -0.6989 0.2056 726.0022 -3.40 0.00071 ***
## normalized_mean_rect -0.0208 0.1700 726.0002 -0.12 0.90257
## normalized_mean_curv -0.4586 0.1523 726.0034 -3.01 0.00270 **
## event1:domain1 -0.0903 0.0868 726.0002 -1.04 0.29868
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.004
## domain1 -0.005 -0.055
## nrmlzd_qr_c -0.087 -0.036 0.084
## nrmlzd_mn_l 0.067 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.056 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.007 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.022 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.058 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.003 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023check_collinearity(APC_model_singlebetas)RFC_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "precentral_inferiorfrontal_R"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
summary(RFC_model_singlebetas)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "precentral_inferiorfrontal_R")
##
## REML criterion at convergence: 3629
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.369 -0.594 0.030 0.630 3.362
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 3.19 1.79
## Residual 5.81 2.41
## Number of obs: 768, groups: subjectID, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.8432 0.3324 32.8136 11.56 4.1e-13 ***
## event1 0.2549 0.0873 727.0000 2.92 0.0036 **
## domain1 0.2856 0.1493 727.0000 1.91 0.0562 .
## normalized_iqr_contrast -0.2156 0.0889 727.0000 -2.43 0.0155 *
## normalized_mean_luminance -0.2577 0.1893 727.0000 -1.36 0.1738
## normalized_mean_motion 0.1398 0.1394 727.0000 1.00 0.3161
## normalized_mean_hsf -0.6089 0.2062 727.0000 -2.95 0.0032 **
## normalized_mean_rect 0.1171 0.1705 727.0000 0.69 0.4926
## normalized_mean_curv -0.3032 0.1528 727.0000 -1.98 0.0475 *
## event1:domain1 -0.0143 0.0871 727.0000 -0.16 0.8696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.005
## domain1 -0.006 -0.055
## nrmlzd_qr_c -0.106 -0.036 0.084
## nrmlzd_mn_l 0.081 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.068 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.009 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.027 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.071 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.003 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023check_collinearity(RFC_model_singlebetas)RSTS_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "superiortemporal_R"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(RSTS_model_singlebetas)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "superiortemporal_R")
##
## REML criterion at convergence: 3402
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.790 -0.573 0.014 0.589 3.773
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 2.50595 1.5830
## extracted_run_number (Intercept) 0.00173 0.0416
## Residual 4.29415 2.0722
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.0826 0.2951 28.1989 7.06 1.1e-07 ***
## event1 0.1045 0.0751 726.1236 1.39 0.16439
## domain1 -0.3892 0.1284 726.0052 -3.03 0.00253 **
## normalized_iqr_contrast -0.1251 0.0765 726.0009 -1.64 0.10224
## normalized_mean_luminance -0.0367 0.1628 726.0027 -0.23 0.82180
## normalized_mean_motion -0.0225 0.1199 726.0149 -0.19 0.85101
## normalized_mean_hsf -0.6320 0.1773 726.0049 -3.56 0.00039 ***
## normalized_mean_rect 0.1491 0.1466 726.0000 1.02 0.30962
## normalized_mean_curv -0.0527 0.1314 726.0078 -0.40 0.68854
## event1:domain1 0.0298 0.0749 725.9999 0.40 0.69066
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.005
## domain1 -0.005 -0.055
## nrmlzd_qr_c -0.103 -0.036 0.084
## nrmlzd_mn_l 0.078 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.066 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.009 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.026 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.069 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.003 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023check_collinearity(RSTS_model_singlebetas)LSTS_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "superiortemporal_L"),
formula = meanbeta ~ event * domain + event_n_within_run +domain * event + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(LSTS_model_singlebetas)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + event_n_within_run + domain * event +
## normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion +
## normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv +
## (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "superiortemporal_L")
##
## REML criterion at convergence: 3719
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.839 -0.511 0.058 0.546 3.911
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.386 2.094
## extracted_run_number (Intercept) 0.111 0.333
## Residual 6.414 2.533
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.1481 0.4877 11.9224 2.35 0.037 *
## event1 0.0625 0.0919 725.0360 0.68 0.496
## domain1 -0.3916 0.1570 724.9993 -2.50 0.013 *
## event_n_within_run -0.0305 0.0102 732.2838 -2.98 0.003 **
## normalized_iqr_contrast -0.0894 0.0934 724.9989 -0.96 0.339
## normalized_mean_luminance 0.1916 0.1990 725.0016 0.96 0.336
## normalized_mean_motion -0.0932 0.1465 725.0002 -0.64 0.525
## normalized_mean_hsf -0.3465 0.2167 725.0001 -1.60 0.110
## normalized_mean_rect -0.2323 0.1792 725.0003 -1.30 0.195
## normalized_mean_curv 0.0402 0.1607 725.0061 0.25 0.803
## event1:domain1 -0.0997 0.0916 725.0082 -1.09 0.276
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 evn___ nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m
## event1 -0.015
## domain1 0.000 -0.055
## evnt_n_wth_ -0.376 0.050 -0.010
## nrmlzd_qr_c -0.072 -0.037 0.084 -0.011
## nrmlzd_mn_l 0.066 -0.027 0.212 -0.021 -0.032
## nrmlzd_mn_m 0.047 0.086 -0.585 0.004 -0.381 -0.212
## nrmlzd_mn_h -0.001 -0.062 0.714 -0.014 0.399 0.344 -0.605
## nrmlzd_mn_r 0.013 -0.009 0.324 0.018 -0.017 -0.432 -0.089
## nrmlzd_mn_c -0.039 -0.056 0.442 -0.031 0.365 0.498 -0.615
## event1:dmn1 0.011 -0.006 0.030 -0.037 0.017 0.002 -0.049
## nrmlzd_mn_h nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## evnt_n_wth_
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r 0.296
## nrmlzd_mn_c 0.497 -0.426
## event1:dmn1 0.030 0.011 0.024check_collinearity(LSTS_model_singlebetas)RSMG_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "supramarginal_R"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(RSMG_model_singlebetas)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "supramarginal_R")
##
## REML criterion at convergence: 3449
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.579 -0.562 0.005 0.565 6.309
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 7.3025 2.702
## extracted_run_number (Intercept) 0.0787 0.281
## Residual 4.3692 2.090
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.5859 0.5250 20.0189 4.93 0.000081 ***
## event1 0.1017 0.0757 726.0185 1.34 0.179
## domain1 0.2934 0.1295 726.0007 2.27 0.024 *
## normalized_iqr_contrast -0.1754 0.0771 726.0001 -2.27 0.023 *
## normalized_mean_luminance -0.4043 0.1642 726.0003 -2.46 0.014 *
## normalized_mean_motion 0.1410 0.1209 726.0021 1.17 0.244
## normalized_mean_hsf -0.4244 0.1788 726.0007 -2.37 0.018 *
## normalized_mean_rect 0.1570 0.1479 725.9999 1.06 0.289
## normalized_mean_curv -0.2804 0.1325 726.0011 -2.12 0.035 *
## event1:domain1 0.1810 0.0755 725.9999 2.40 0.017 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.003
## domain1 -0.003 -0.055
## nrmlzd_qr_c -0.058 -0.036 0.084
## nrmlzd_mn_l 0.044 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.037 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.005 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.015 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.039 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.002 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023plot(allEffects(RSMG_model_singlebetas))
lsmeans(RSMG_model_singlebetas, pairwise ~ event | domain)$contrasts## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.565 0.213 726 2.649 0.0080
##
## domain = psychology-action:
## contrast estimate SE df t.ratio p.value
## unexp - exp -0.158 0.214 726 -0.739 0.4600
##
## Degrees-of-freedom method: kenward-rogercheck_collinearity(RSMG_model_singlebetas)LSMG_model_singlebetas <- lmer(data = single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final == "supramarginal_L"),
formula = meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|extracted_run_number) + (1|subjectID))
summary(LSMG_model_singlebetas)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## meanbeta ~ event * domain + normalized_iqr_contrast + normalized_mean_luminance +
## normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect +
## normalized_mean_curv + (1 | extracted_run_number) + (1 | subjectID)
## Data: single_betas_focal_psychphys_runs12 %>% filter(ROI_name_final ==
## "supramarginal_L")
##
## REML criterion at convergence: 3379
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.060 -0.600 0.012 0.579 5.640
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 4.11543 2.0287
## extracted_run_number (Intercept) 0.00379 0.0616
## Residual 4.07805 2.0194
## Number of obs: 768, groups: subjectID, 32; extracted_run_number, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.4037 0.3715 30.2024 6.47 3.7e-07 ***
## event1 0.0665 0.0731 726.1058 0.91 0.3636
## domain1 0.2683 0.1252 726.0045 2.14 0.0324 *
## normalized_iqr_contrast -0.0619 0.0745 726.0008 -0.83 0.4064
## normalized_mean_luminance -0.4990 0.1587 726.0024 -3.14 0.0017 **
## normalized_mean_motion -0.0931 0.1168 726.0128 -0.80 0.4256
## normalized_mean_hsf -0.0259 0.1728 726.0042 -0.15 0.8808
## normalized_mean_rect 0.2933 0.1429 726.0001 2.05 0.0405 *
## normalized_mean_curv -0.0639 0.1280 726.0067 -0.50 0.6176
## event1:domain1 0.2042 0.0730 726.0000 2.80 0.0053 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 nrmlzd_q_ nrmlzd_mn_l nrmlzd_mn_m nrmlzd_mn_h
## event1 0.004
## domain1 -0.004 -0.055
## nrmlzd_qr_c -0.079 -0.036 0.084
## nrmlzd_mn_l 0.061 -0.026 0.212 -0.033
## nrmlzd_mn_m 0.051 0.086 -0.585 -0.381 -0.212
## nrmlzd_mn_h -0.007 -0.062 0.714 0.399 0.344 -0.605
## nrmlzd_mn_r 0.020 -0.010 0.325 -0.017 -0.432 -0.089 0.296
## nrmlzd_mn_c -0.053 -0.054 0.442 0.365 0.497 -0.615 0.497
## event1:dmn1 -0.003 -0.004 0.029 0.017 0.002 -0.049 0.029
## nrmlzd_mn_r nrmlzd_mn_c
## event1
## domain1
## nrmlzd_qr_c
## nrmlzd_mn_l
## nrmlzd_mn_m
## nrmlzd_mn_h
## nrmlzd_mn_r
## nrmlzd_mn_c -0.425
## event1:dmn1 0.012 0.023plot(allEffects(LSMG_model_singlebetas))
lsmeans(LSMG_model_singlebetas, pairwise ~ event | domain)$contrasts## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.541 0.206 726 2.626 0.0088
##
## domain = psychology-action:
## contrast estimate SE df t.ratio p.value
## unexp - exp -0.275 0.207 726 -1.330 0.1838
##
## Degrees-of-freedom method: kenward-rogercheck_collinearity(LSMG_model_singlebetas)7.5 PART 5b: Per task event analysis
single_betas_pertask_data <- single_betas %>%
filter(event != "fam") %>%
filter(str_length(extracted_run_number) == 4)
tasks <- unique(single_betas_pertask_data$task)
ROIs <- unique(single_betas_pertask_data$ROI_name_final)
single_betas_pertask_data$event <- relevel(single_betas_pertask_data$event, ref = "unexp")
runs <- c("all_runs", "runs12")single_betas_pertask_results <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (task in tasks) {
cur_task <- task
for (run in runs) {
if (run == "all_runs") {
ROIdata <-
single_betas_pertask_data %>% filter(ROI_name_final == curROI,
task == cur_task)
}
else if (run == "runs12") {
ROIdata <-
single_betas_pertask_data %>%
filter(ROI_name_final == curROI,
task == cur_task) %>%
filter(extracted_run_number == "run1" |
extracted_run_number == "run2")
}
# print("Current ROI: ", curROI)
model <- lmer(data = ROIdata,
formula = meanbeta ~ event + (1 | subjectID)) # had to remove run random int (convergence))
focal_region <- ROIdata$focal_region[1]
task_domain <- as.character(ROIdata$domain[1])
ROI_category <- ROIdata$ROI_category[1]
model_null <-
update(model, formula = ~ . - event) # Without interaction term
BF_BIC_event <-
data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_null) - BIC(model)
) / 2))) # BICs to Bayes factor
# BF_initial <- rbind(BF_BIC_event)
rownames(BF_BIC_event) <- c("event1")
BFs <- rbind(BF_intercept_row,
BF_BIC_event)
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <-
rbind(effect_sizes_intercept_row, effect_sizes_initial)
cur_top_voxels_N <- 100
curRun <- run
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
focal_region,
curROI,
cur_task,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:4, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"task",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d",
"BF",
"BF_interpretation"
)
single_betas_pertask_results <-
rbind(single_betas_pertask_results, modelsummary)
# model_parameters(model) %>% print_md()
cat(
tab_model(
model,
show.std = TRUE,
show.se = TRUE,
show.stat = TRUE,
auto.label = TRUE,
title = paste0(
"Exploratory Univariate ROI results for the task ",
cur_task,
"in the region: ",
curROI,
"in runs ",
run
)
)$knitr,
" \n \n"
)
}
}
}single_betas_pertask_results <- single_betas_pertask_results %>%
mutate_at(c(8:16), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
single_betas_pertask_results$effect <-
as.factor(single_betas_pertask_results$effect)
levels(single_betas_pertask_results$effect) <-
c("intercept", "event")
rownames(single_betas_pertask_results) <- NULL
# write_rds(single_betas_pertask_results, here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds"))single_betas_pertask_results <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds"))
DT::datatable(dplyr::arrange(single_betas_pertask_results, extracted_run_number, effect), options = list(scrollX = TRUE))8 Part 6: Overlap of top 100 voxels
overlap_data_prelim <- read_csv(here("input_data/ROI_overlap_outputs_Apr19_2023.csv")) %>%
rename(ROI_2 = roi_2,
subjID = participantID)## Parsed with column specification:
## cols(
## participantID = col_character(),
## ROI_1 = col_character(),
## roi_2 = col_character(),
## overlap_dices_coef = col_double()
## )MD_physics_APC_bilateral <- overlap_data_prelim %>%
filter((
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" &
ROI_2 == "MD_ROI_3_13_bilateral"
) |
(
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" &
ROI_2 == "MD_ROI_3_13_bilateral"
)
) %>%
mutate(hemisphere = "bilateral")
MD_physics_APC_byhem <- overlap_data_prelim %>%
filter((
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" &
ROI_2 == "MD_ROI_13"
) |
(
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" &
ROI_2 == "MD_ROI_3"
)
) %>%
mutate(hemisphere = case_when(ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-right_inflated_notSTS" ~ "right",
ROI_1 == "DOTS_VOE234_phys_p20_cluster-4-left_inflated_notSTS" ~ "left"))
names(MD_physics_APC_bilateral)## [1] "subjID" "ROI_1" "ROI_2"
## [4] "overlap_dices_coef" "hemisphere"names(MD_physics_APC_byhem)## [1] "subjID" "ROI_1" "ROI_2"
## [4] "overlap_dices_coef" "hemisphere"MD_physics_APC_bilateral <- rbind(MD_physics_APC_bilateral, MD_physics_APC_byhem)overlap_summary <- ggplot(MD_physics_APC_bilateral, aes(hemisphere, overlap_dices_coef)) +
geom_boxplot() +
geom_point(fun = "mean", stat = "summary", shape = 18, colour = "red", size = 4) +
geom_jitter(height = 0, width = .1) +
ylab("Dice's Coefficient") +
xlab("APC ROI (Bilateral = Confirmatory)")
# geom_text_repel(aes(label = subjID), size = 3) +
# ggtitle("Dice's Coefficient between LSMG and RSMG \n to bilateral APC, left APC, and right APC")overlap_summary
APC_SMG_overlap_results <- MD_physics_APC_bilateral %>%
group_by(hemisphere) %>%
summarise(mean_DC = mean(overlap_dices_coef),
median_DC = median(overlap_dices_coef),
range_low = range(overlap_dices_coef)[1],
range_high = range(overlap_dices_coef)[2])
MD_physics_APC_bilateral %>%
filter(overlap_dices_coef == 0) %>%
group_by(hemisphere) %>%
summarise(n = n())MD_physics_APC_bilateral %>%
filter(overlap_dices_coef != 0,
hemisphere == "bilateral") %>%
group_by(subjID) MD_physics_APC_byhem %>%
group_by(hemisphere) %>%
summarise(mean_DC = mean(overlap_dices_coef),
median_DC = median(overlap_dices_coef),
range_low = range(overlap_dices_coef)[1],
range_high = range(overlap_dices_coef)[2])MD_physics_APC_byhem %>%
filter(overlap_dices_coef == 0) %>%
group_by(hemisphere) %>%
summarise(n = n())# write_rds(APC_SMG_overlap_results, path = here("outputs/univariate_results/APC_SMG_overlap_results.Rds"))9 Part 6: Effects per task
single_betas_pertask_results <- read_rds(here("outputs/univariate_results/univar_results_top100_singlebetas_pertask.Rds"))
# filter((domain == "both" & extracted_run_number == "all_runs") | (domain != "both" & extracted_run_number == "runs12"))tasks <- c("infer-constraint", "agent-solidity")
single_betas_pertask_results_focal <- single_betas_pertask_results %>%
filter(focal_region == 1,
effect == "event") %>%
filter((extracted_run_number == "all_runs" & task %in% tasks) | (extracted_run_number == "runs12" & !task %in% tasks)) %>%
select(extracted_run_number, ROI, task, domain, d, BF, B, SE, CI_95_lower, CI_95_upper, star)
single_betas_pertask_results_focal$task <- factor(single_betas_pertask_results_focal$task, levels = c("solidity", "permanence", "goal", "efficiency", "agent-solidity", "infer-constraint"))
single_betas_pertask_results_focal$ROI <- factor(single_betas_pertask_results_focal$ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_bilateral", "MT_bilateral"))
pertask_b_pertask <- ggplot(single_betas_pertask_results_focal,
aes(ROI, B, colour = task, group = task, label = star)) +
geom_point() +
geom_errorbar(
aes(ymin = B - SE, ymax = B + SE,
colour = task),
position = position_dodge(width = .9),
width = .2
) +
geom_text(size = 4, aes(ROI, B + SE + .02)) +
# geom_line() +
geom_hline(yintercept = 0) +
facet_wrap( ~ task, nrow = 1) +
scale_colour_manual(values = c("#25a88b", "#19626a", "#f24b00", "#deae1e", "#6600cd", "#fba0d0")) +
ylab("B and SE") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
coord_cartesian(ylim = c(-.6, .8))
pertask_b_perROI <- ggplot(single_betas_pertask_results_focal,
aes(task, B, colour = task, group = task, label = star)) +
geom_point() +
geom_errorbar(
aes(ymin = B - SE, ymax = B + SE,
colour = task),
position = position_dodge(width = .9),
width = .2
) +
geom_text(size = 4, aes(task, B + SE + .02)) +
# geom_line() +
geom_hline(yintercept = 0) +
facet_wrap( ~ ROI, nrow = 1) +
scale_colour_manual(values = c("#25a88b", "#19626a", "#f24b00", "#deae1e", "#6600cd", "#fba0d0")) +
ylab("B and SE") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
coord_cartesian(ylim = c(-.6, .8))pertask_b_pertask
pertask_b_perROI
10 Part 7: Region by region univar effect size
This analysis accounts for low level visual confounds.
data_for_modeling <- single_betas %>%
filter(event != "fam",
top_n_voxels == "100",
domain != "psychology-environment") %>%
filter(extracted_run_number == "run1" | extracted_run_number == "run2") %>%
filter(!str_detect(ROI_name_final, "bilateral"))
data_for_modeling$ROI_name_final <- as.factor(data_for_modeling$ROI_name_final)
ROIs <- sort(unique(data_for_modeling$ROI_name_final))
data_for_modeling$event <- relevel(data_for_modeling$event, ref = "unexp")modelsummaries_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
# FIRST COMPUTE EVENT EFFECT SIZES, FOR EACH DOMAIN
physics_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
domain == "physics")
psychology_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
domain == "psychology-action")
# had to remove random intercept for runs (convergence))
physics_model <- lmer(data = physics_data,
formula = meanbeta ~ event + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
psychology_model <- lmer(data = psychology_data,
formula = meanbeta ~ event + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
physics_event_effect_size <- lme.dscore(physics_model, physics_data, "lme4") %>%
select(d) %>%
slice(1)
psychology_event_effect_size <- lme.dscore(psychology_model, psychology_data, "lme4") %>%
select(d) %>%
slice(1)
# SECOND COMPUTE DOMAIN EFFECT SIZES, FOR EACH EVENT
expected_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
event == "exp")
unexpected_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
event == "unexp")
# had to remove random intercept for runs (convergence))
exp_model <- lmer(data = expected_data,
formula = meanbeta ~ domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
unexp_model <- lmer(data = unexpected_data,
formula = meanbeta ~ domain + normalized_iqr_contrast + normalized_mean_luminance + normalized_mean_motion + normalized_mean_hsf + normalized_mean_rect + normalized_mean_curv + (1|subjectID))
exp_domain_effect_size <-
lme.dscore(exp_model, expected_data, "lme4") %>%
select(d) %>%
slice(1)
unexp_domain_effect_size <-
lme.dscore(unexp_model, unexpected_data, "lme4") %>%
select(d) %>%
slice(1)
if (physics_data$ROI_category[1] == "MD" | physics_data$ROI_category[1] == "early_visual") {
ROI_domain <- "general"
} else {
ROI_domain <- "specific"
}
modelsummary <-
cbind(
ROI_domain,
physics_data$ROI_category[1],
curROI,
physics_event_effect_size,
psychology_event_effect_size,
exp_domain_effect_size,
unexp_domain_effect_size
) %>%
as.data.frame()
modelsummaries_init <-
rbind(modelsummaries_init, modelsummary)
}colnames(modelsummaries_init) <-
c("ROI_domain",
"ROI_category",
"ROI",
"physics_event_effect_size",
"psychology_event_effect_size",
"exp_domain_effect_size",
"unexp_domain_effect_size"
)
rownames(modelsummaries_init) <- NULL
# write_rds(modelsummaries_init, here("outputs/univariate_results/univar_manyregions_viscontrols_results.Rds"))univar_manyregions_results <- read_rds( here("outputs/univariate_results/univar_manyregions_viscontrols_results.Rds"))plot1 <- ggplot(univar_manyregions_results, aes(psychology_event_effect_size, physics_event_effect_size)) +
geom_hline(yintercept = 0, linetype = 'dotted') +
geom_vline(xintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-0.5, 1), ylim = c(-0.5, 1)) +
geom_point(aes(colour = ROI_category), size = 3) +
geom_smooth(method = "lm") +
facet_wrap(~ROI_domain, scales = "free") +
scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
# coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
xlab("VOE Effect Size (d), Psychology Events") +
ylab("VOE Effect Size (d), Physics Events") +
ggtitle(label = "Exp 2 VOE effects for each domain across regions") +
geom_text_repel(aes(label = ROI), size = 3) +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.position = "none")
plot2 <- ggplot(univar_manyregions_results, aes(exp_domain_effect_size, unexp_domain_effect_size)) +
geom_hline(yintercept = 0, linetype = 'dotted') +
geom_vline(xintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-1.2, 1.2), ylim = c(-1.2, 1.2)) +
geom_point(aes(colour = ROI_category), size = 3) +
geom_smooth(method = "lm") +
facet_wrap(~ROI_domain, scales = "free") +
scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
# coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
xlab("Domain Effect Size (d), Expected Events") +
ylab("Domain Effect Size (d), Unexpected Events") +
ggtitle(label = "Exp 2 Domain effects for each event across regions") +
geom_text_repel(aes(label = ROI), size = 3) +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.position = "none")plot1## `geom_smooth()` using formula = 'y ~ x'
plot2## `geom_smooth()` using formula = 'y ~ x'
# instead of testing whether the linear relationship between x and y is 0,
# we test for independence instead. H0 is that x and y are independent; F_{XY}(x,y) = F_X(x) F_Y(y).
observed_cors_ind <- univar_manyregions_results %>%
group_by(ROI_domain) %>%
summarise(
cor_domain_within_event =
np.cor.test(
exp_domain_effect_size,
unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE
)$estimate,
p_domain_within_event =
np.cor.test(
exp_domain_effect_size,
unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE
)$p.value,
cor_event_within_domain =
np.cor.test(
psychology_event_effect_size,
physics_event_effect_size,
alternative = "two.sided",
independent = TRUE
)$estimate,
p_event_within_domain =
np.cor.test(
psychology_event_effect_size,
physics_event_effect_size,
alternative = "two.sided",
independent = TRUE
)$p.value
)observed_cors_ind# write_rds(observed_cors_ind, path = here("outputs/univariate_results/study2_univar_manyregions_observed_cor.Rds"))10.0.1 Bootstrap stats
set.seed(2020)
## First define a function to work out the difference of corrs:
diff_corr <- function(data, indices) {
data <- data[indices,]
cor1 <-
np.cor.test(DS_results$exp_domain_effect_size,
DS_results$unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE,
parallel = TRUE)$estimate
cor2 <-
np.cor.test(
data$psychology_event_effect_size,
data$physics_event_effect_size,
alternative = "two.sided",
independent = TRUE,
parallel = TRUE)$estimate
return(cor1 - cor2)
}10.0.2 Domain specific
DS_results <- univar_manyregions_results %>%
filter(ROI_domain == "specific")doMC::registerDoMC(cores = 2) # for running in parallel
# Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DS <- boot(
# data = DS_results,
# R = 4000,
# statistic = diff_corr,
# stype = "i"
# )
# saveRDS(res_boot_DS, here("outputs/univariate_results/study2_univariate_dsregions_4000perms_confirmatory.Rds"))
res_boot_DS <- read_rds(here("outputs/univariate_results/study2_univariate_dsregions_4000perms_confirmatory.Rds"))
## Retrieve the empirical 95% confidence interval:
ds_results <- append(boot.ci(res_boot_DS, type = "perc", conf = 0.95),boot.pval(res_boot_DS))
# saveRDS(ds_results, here("outputs/univariate_results/study2_univariate_dsregions_summary.Rds"))10.0.3 Domain general
DG_results <- univar_manyregions_results %>%
filter(ROI_domain == "general")## Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DG <- boot(data = DG_results,
# R = 4000,
# statistic = diff_corr,
# stype = "i")
# saveRDS(res_boot_DG, here("outputs/univariate_results/study2_univariate_dgregions_4000perms_confirmatory.Rds"))
res_boot_DG <- readRDS(here("outputs/univariate_results/study2_univariate_dgregions_4000perms_confirmatory.Rds"))
plot(res_boot_DG)
## Retrieve the empirical 95% confidence interval (adjusted to 90% because of one-tailed prediction):
dg_results <- append(boot.ci(res_boot_DG, type = "perc", conf = 0.95),boot.pval(res_boot_DG))
# saveRDS(dg_results, here("outputs/univariate_results/study2_univariate_dgregions_summary.Rds"))